# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
# install.packages(c("tidyverse", "estimatr", "coefplot", "reshape2))

library(tidyverse)
library(estimatr)
library(coefplot)
library(reshape2)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")

source("opeds_source.R")

dvs <-
  c(
    "dv_amtrak_main_w1",
    "dv_climate_main_w1",
    "dv_flat_main_w1",
    "dv_vets_main_w1",
    "dv_wall_main_w1",
    "dv_amtrak_scale_w1",
    "dv_climate_scale_w1",
    "dv_flat_scale_w1",
    "dv_vets_scale_w1",
    "dv_wall_scale_w1"
  )

dv_labels <-
  rep(c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street"), 2)

Zs <-
  rep(c("amtrak", "climate", "paul", "veterans", "wallstreet"), 2)

parties <- c("Democrat", "Independent", "Republican")

df <- NULL

for (j in 1:length(Zs)) {
  for (k in 1:length(parties)) {
    local_subset <-
      mturk_opeds$Z %in% c("control", Zs[j]) &
      mturk_opeds$pid_3_cat == parties[k]
    local_df <- mturk_opeds[local_subset, ]
    local_formula <- formula(paste0(dvs[j], "~Z"))
    
    summary_df <-
      lm_robust(local_formula, data = local_df) %>%
      tidy %>%
      mutate(
        treatment = Zs[j],
        dependent_variable = dv_labels[j],
        dv = dvs[j],
        party = parties[k],
        sample = "Mechanical Turk"
      )
    
    df <- bind_rows(df, summary_df)
  }
}


dvs <- c("dv_amtrak_main_w1", "dv_flat_main_w1", "dv_vets_main_w1", "dv_wall_main_w1",
         "dv_amtrak_scale_w1", "dv_flat_scale_w1", "dv_vets_scale_w1", "dv_wall_scale_w1")
dv_labels <- rep(c("Amtrak", "Flat Tax", "Veterans", "Wall Street"), 2)
Zs <- rep(c("amtrak", "paul", "veterans", "wallstreet"), 2)
parties <- c("Democrat", "Independent", "Republican")

for(j in 1:length(Zs)) {
  for (k in 1:length(parties)) {
    local_subset <-
      elite_opeds$Z %in% c("control", Zs[j]) &
      elite_opeds$pid_3_cat == parties[k]
    local_df <- elite_opeds[local_subset, ]
    local_formula <- formula(paste0(dvs[j], "~Z"))

    
    summary_df <-
      lm_robust(local_formula, data = local_df) %>%
      tidy %>%
      mutate(
        treatment = Zs[j],
        dependent_variable = dv_labels[j],
        dv = dvs[j],
        party = parties[k],
        sample = "Elites"
      )
    
    df <- bind_rows(df, summary_df)
  }
}


df <-
  df %>%
  filter(term != "(Intercept)") %>%
  mutate(dv_type = ifelse(grepl("scale", dv), "Scale DV", "Main DV"))


g <- 
  ggplot(df, aes(y = dependent_variable, x = estimate, group = party, color = party, shape = party)) + 
  geom_point(position = position_dodgev(height = .5), size = 2) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 position = position_dodgev(height = .5),
                 height = 0) +
  scale_color_manual(values = c("blue", "purple", "red")) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    axis.title = element_blank(),
    legend.key.width = unit(4, "lines"),
    strip.background = element_blank()
  ) +
  ylab("Estimated Treatment Effects") +
  coord_cartesian(xlim = c(-.1, 1.5)) +
  #ggtitle("Effects of Treatment, by Party and Experimental Sample") +
  facet_grid(sample ~ dv_type, scales = "free") 

g




# Precision weighted averages ---------------------------------------------

df %>%
  group_by(dv_type, sample, party) %>%
  summarise(average_ate = weighted.mean(estimate, 1/(std.error^2)),
            average_ate_se = sqrt(1/sum(1/std.error^2)))

# ftests ------------------------------------------------------------------

dvs <-
  c(
    "dv_amtrak_main_w1",
    "dv_climate_main_w1",
    "dv_flat_main_w1",
    "dv_vets_main_w1",
    "dv_wall_main_w1",
    "dv_amtrak_scale_w1",
    "dv_climate_scale_w1",
    "dv_flat_scale_w1",
    "dv_vets_scale_w1",
    "dv_wall_scale_w1"
  )
dv_labels <-
  rep(c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street"), 2)
Zs <-
  rep(c("amtrak", "climate", "paul", "veterans", "wallstreet"), 2)

df <- NULL

for(j in 1:length(Zs)){
  local_subset <- mturk_opeds$Z %in% c("control",Zs[j])
  local_df <- mturk_opeds[local_subset,]
  local_formula_r <- formula(paste0(dvs[j], "~Z + pid_3_cat"))
  local_formula_u <- formula(paste0(dvs[j], "~Z*pid_3_cat"))
  fit_r <- lm_robust(local_formula_r, data = local_df)
  fit_u <- lm_robust(local_formula_u, data = local_df)
  anova_fit <- waldtest(fit_r, fit_u, test = "F")
  
  anova_fit$`Pr(>F)`[2]
  
  summary_df <- data.frame(
    treatment = Zs[j],
    dependent_variable = dv_labels[j],
    dv = dvs[j],
    p = anova_fit$`Pr(>F)`[2],
    sample = "Mechanical Turk"
  )
  df <- bind_rows(df, summary_df)
}

# elites ------------------------------------------------------------------

dvs <- c("dv_amtrak_main_w1", "dv_flat_main_w1", "dv_vets_main_w1", "dv_wall_main_w1",
         "dv_amtrak_scale_w1", "dv_flat_scale_w1", "dv_vets_scale_w1", "dv_wall_scale_w1")
dv_labels <- rep(c("Amtrak", "Flat Tax", "Veterans", "Wall Street"), 2)
Zs <- rep(c("amtrak", "paul", "veterans", "wallstreet"), 2)

for(j in 1:length(Zs)){
  local_subset <- elite_opeds$Z %in% c("control", Zs[j])
  local_df <- elite_opeds[local_subset, ]
  local_formula_r <- formula(paste0(dvs[j], "~Z + pid_3_cat"))
  local_formula_u <- formula(paste0(dvs[j], "~Z*pid_3_cat"))
  fit_r <- lm_robust(local_formula_r, data = local_df)
  fit_u <- lm_robust(local_formula_u, data = local_df)
  
  anova_fit <- waldtest(fit_r, fit_u, test = "F")
  
  anova_fit$`Pr(>F)`[2]
  
  summary_df <- data.frame(
    treatment = Zs[j],
    dependent_variable = dv_labels[j],
    dv = dvs[j],
    p = anova_fit$`Pr(>F)`[2],
    sample = "Elites"
  )
  df <- bind_rows(df, summary_df)
}


mturk_ftests <-
  df %>%
  mutate(dv_type = ifelse(grepl("main", dv), "Main DV", "Scale DV")) %>%
  filter(sample == "Mechanical Turk") %>%
  dcast(dependent_variable ~ dv_type, value.var = "p")

mturk_ftests

elite_ftests <- 
  df %>% 
  mutate(dv_type = ifelse(grepl("main", dv), "Main DV", "Scale DV")) %>%
  filter(sample == "Elites") %>%
  dcast(dependent_variable ~ dv_type, value.var = "p")

elite_ftests
